Estimate Treatment Effects RCT

Agustin Calatroni
AXC

Tuesday, the 31^th of January, 2023

Outline

How to estimate treatment effects in longitudinal RCT’s
Problem
Data
Models R &SAS
Results
How to extend these estimation
How to aggregate & communicate these results

Different ways to estimate treatment effects in randomised controlled trials

Abstract

Regarding the analysis of longitudinal RCT data there is a debate going on whether an adjustment for the baseline value of the outcome variable should be made. When an adjustment is made, there is a lot of misunderstanding regarding the way this should be done. Therefore, the aims to:

to explain different methods used to estimate treatment effects in RCTs
to illustrate the different methods with a real data example
to give an advise on how to analyse RCT data

ADaM Basic Data Structure (BDS)

Example dataset

Treatment of Lead Exposed Children Trial (N=100)

Wide
Long

id	trt	w_0	w_1	w_4	w_6
1	Placebo	31.9	26.9	25.8	23.8
2	Active	26.5	14.8	19.5	21.0
3	Active	25.8	23.0	19.1	23.2
4	Placebo	26.5	24.5	22.0	22.5
5	Active	20.4	2.8	3.2	9.4

99	Active	21.9	7.6	10.8	13.0
100	Active	20.7	8.1	25.7	12.3

id	trt	week	lead
1	Placebo	0	31.9
1	Placebo	1	26.9
1	Placebo	4	25.8
1	Placebo	6	23.8
2	Active	0	26.5

100	Active	4	25.7
100	Active	6	12.3

Table (wo/missing)

Treatment of Lead Exposed Children Trial (N=100)

Code

tbl_summary(data = tlc_w,
            by = trt, 
            include = starts_with('w'),
            digit = list(everything() ~ c(1,2)),
            statistic = list(all_continuous() ~ "{mean} ({sd})"),
            label = list( w_0 ~ 'Week 0',
                          w_1 ~ 'Week 1',
                          w_4 ~ 'Week 4',
                          w_6 ~ 'Week 6')) %>% 
   add_difference(test = list(all_continuous() ~ "t.test"),
                  pvalue_fun = ~style_pvalue(.x, digits = 2),
                  estimate_fun = list(all_continuous() ~ function(x) style_number(x, digits=1)) ) %>% 
   bold_labels() %>% 
   modify_column_merge(pattern = "{estimate} ({ci})") %>% 
   modify_header(label = '',
                 estimate = '**Difference (95% CI)**',
                 all_stat_cols(FALSE) ~ "**{level}**<br>n = {n}" ) %>% 
   modify_spanning_header(all_stat_cols(FALSE) ~ "**Treatment**") %>% 
   # modify_footnote(everything() ~ NA) %>%
   add_n()

	N	Treatment		Difference (95% CI)²	p-value²
	N	Active n = 50¹	Placebo n = 50¹	Difference (95% CI)²	p-value²
Week 0	100	26.5 (5.02)	28.0 (4.98)	-1.5 (-3.4, 0.5)	0.15
Week 1	100	13.5 (7.67)	24.7 (5.46)	-11.1 (-13.8, -8.5)	<0.001
Week 4	100	15.5 (7.85)	24.1 (5.75)	-8.6 (-11.3, -5.8)	<0.001
Week 6	100	20.8 (9.25)	23.6 (5.64)	-2.9 (-5.9, 0.2)	0.063
¹ Mean (SD)
² Welch Two Sample t-test

Tables (w/missing)

Treatment of Lead Exposed Children Trial with Missing Data (N=100)

TABLE-1
TABLE-2

Code

tbl_summary(data = tlcmiss_w ,
            by = trt, 
            missing = 'no',
            include = starts_with('w'),
            digit = list(everything() ~ c(1,2)),
            statistic = list(all_continuous() ~ "{mean} ({sd})"),
            label = list( w_0 ~ 'Week 0',
                          w_1 ~ 'Week 1',
                          w_4 ~ 'Week 4',
                          w_6 ~ 'Week 6')) %>% 
   add_difference(test = list(all_continuous() ~ "t.test"),
                  pvalue_fun = ~style_pvalue(.x, digits = 2),
                  estimate_fun = list(all_continuous() ~ function(x) style_number(x, digits=1)) ) %>% 
   bold_labels() %>% 
   modify_column_merge(pattern = "{estimate} ({ci})") %>% 
   modify_header(label = '',
                 estimate = '**Difference (95% CI)**',
                 all_stat_cols(FALSE) ~ "**{level}**<br>n = {n}" ) %>% 
   modify_spanning_header(all_stat_cols(FALSE) ~ "**Treatment**") %>% 
   # modify_footnote(everything() ~ NA) %>%
   add_n()

	N	Treatment		Difference (95% CI)²	p-value²
	N	Active n = 50¹	Placebo n = 50¹	Difference (95% CI)²	p-value²
Week 0	100	26.5 (5.02)	28.0 (4.98)	-1.5 (-3.4, 0.5)	0.15
Week 1	88	13.4 (8.11)	24.5 (5.46)	-11.1 (-14.1, -8.1)	<0.001
Week 4	81	15.2 (7.08)	24.1 (5.80)	-8.9 (-11.8, -5.9)	<0.001
Week 6	74	19.3 (10.48)	23.7 (5.69)	-4.4 (-8.9, 0.1)	0.055
¹ Mean (SD)
² Welch Two Sample t-test

Code

tbl_reg_summary(data = tlcmiss_w,
                by = trt, 
                include = starts_with('w'),
                digits = list(everything() ~ c(0,1,2,1,1,1,1,1,0))) %>% 
   add_overall() %>% 
   bold_labels() %>% 
   modify_header(all_stat_cols(TRUE) ~ "**{level}**" )

Characteristic	Overall	Active	Placebo
Week 0
N	100	50	50
Mean (SD)	27.3 (5.03)	26.5 (5.02)	28.0 (4.98)
Median (IQR)	26.6 (23.0, 30.7)	26.2 (22.1, 29.6)	27.1 (23.6, 31.8)
Range	19.7, 41.1	19.7, 41.1	20.7, 40.7
N missing	0	0	0
Week 1
N	88	41	47
Mean (SD)	19.3 (8.78)	13.4 (8.11)	24.5 (5.46)
Median (IQR)	20.9 (12.9, 25.2)	11.6 (6.5, 17.5)	23.9 (20.9, 27.8)
Range	2.8, 40.8	2.8, 39.0	14.9, 40.8
N missing	12	9	3
Week 4
N	81	35	46
Mean (SD)	20.2 (7.73)	15.2 (7.08)	24.1 (5.80)
Median (IQR)	21.0 (15.9, 25.1)	15.1 (9.1, 20.4)	22.5 (19.8, 27.4)
Range	3.0, 38.6	3.0, 27.0	15.3, 38.6
N missing	19	15	4
Week 6
N	74	26	48
Mean (SD)	22.2 (7.93)	19.3 (10.48)	23.7 (5.69)
Median (IQR)	21.2 (18.4, 24.8)	18.5 (13.5, 22.1)	22.4 (20.1, 27.6)
Range	4.1, 63.9	4.1, 63.9	13.5, 43.3
N missing	26	24	2

Code

ggplot(data = tlc_l %>% 
          mutate(week = as.character(week) %>% as.numeric()),
       aes(x = week, y = lead, color = trt)) +
   geom_quasirandom(dodge.width=1, alpha = 0.5, size = 2.5) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   scale_y_continuous(name = 'Lead (micrograms/dL)',
                      limits = c(0, 65)) +
   scale_x_continuous(name = 'Week',
                      breaks = c(0, 1, 4, 6),
                      labels = c('0\nRandomization', '1', '4', '6')) +
   theme_bw(base_size = 20) + 
   theme(panel.grid.minor = element_blank(),
         legend.position = c(0.02, 0.98),
         legend.justification = c('left','top'))

Code

m2 <- mmrm(lead ~ trt * week + us(week | id), data = tlc_l)
e2 <- emmeans(m2, pairwise ~ trt | week) %>% 
   as.data.frame() %>% 
   mutate(week = as.numeric(week))

f_a <- ggplot(data = tlc_l %>% 
                 mutate(week = as.character(week) %>% as.numeric()),
              aes(x = week, y = lead, color = trt) ) +
   geom_quasirandom(dodge.width=1, alpha = 0.2, size = 2.5, show.legend = FALSE) +
   geom_boxplot(aes(group = interaction(week, trt)),
                fill = 'transparent',
                outlier.shape = NA,
                show.legend = FALSE) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   scale_y_continuous(name = 'Lead (micrograms/dL)',
                      limits = c(0, 65)) +
   scale_x_continuous(name = NULL,
                      breaks = c(0, 1, 4, 6),
                      labels = c('0\nRandomization', '1', '4', '6')) +
   stat_summary(fun = mean,
                geom = "point",
                size = 3,
                position = position_dodge(width = 1)) +
   stat_summary(fun = mean,
                geom = "line",
                size = 1,
                position = position_dodge(width = 1),
                show.legend = FALSE) +
   theme_bw(base_size = 20) +
   theme(panel.grid.minor = element_blank(),
         legend.position = c(0.02, 0.98),
         legend.justification = c('left','top'))

f_b <- ggplot(data = e2 %>% filter(contrast != '.'),
              aes(x = week, y = emmean, group = contrast)) + 
   
   geom_hline(yintercept  = 0, color = 'gray50') +
   geom_point(size = 3) +
   geom_line(size = 1) +
   geom_linerange(aes(ymin = lower.CL,
                      ymax = upper.CL),
                  size = 1) +
   geom_text(aes(y = 4.75,
                 label = gtsummary::style_pvalue(p.value, digits = 2, prepend_p = TRUE)),
             size = 4) +
   scale_color_discrete(name = 'Treatment') +
   scale_y_continuous(name = 'Difference',
                      limits = c(-15, 5)) +
   scale_x_continuous(name = 'Week',
                      expand = c(0.12,0.12),
                      breaks = c(0, 1, 4, 6),
                      labels = c('0\nRandomization', '1', '4', '6')) +
   theme_bw(base_size = 20) +
   theme(panel.grid.minor = element_blank())

f_a / f_b + plot_layout( heights = c(3, 1))

Code

ggplot(data = tlc_l %>% 
          mutate(week = as.character(week) %>% as.numeric()),
       aes(x = week, y = lead, color = trt, group = id)) +
   geom_quasirandom(dodge.width = 0.5, alpha = 0.5, size = 2.5) +
   geom_line(alpha = 0.5,
             position=position_quasirandom(dodge.width = 0.5)) +
   facet_wrap(. ~ trt ) +
   scale_colour_brewer(palette = "Set1",
                       name = 'Treatment') +
   scale_y_continuous(name = 'Lead (micrograms/dL)',
                      limits = c(0, 65)) +
   scale_x_continuous(name = 'Week',
                      breaks = c(0, 1, 4, 6),
                      labels = c('0\nRandomization', '1', '4', '6')) +
   theme_bw(base_size = 18) +
   theme(panel.grid.minor = element_blank(),
         legend.position="none")

trelliscopejs: scalable, flexible & interactive approach to visualizing data

Scatterplot Matrix

Mixed models for repeated measures (MMRM)

Method 1: Longitudinal analysis of covariance

Overall TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0}\]

Code

tlc3_l <- tlc3_l %>% mutate(trt = fct_rev(trt))

m1a <- mmrm(lead ~ trt + base + 
               us(week|id),     
            data = tlc3_l %>% filter(week != '0'))

e1a <- emmeans(m1a, revpairwise ~ trt) %>% 
   as.data.frame() 

summary(m1a)

mmrm fit

Formula:     lead ~ trt + base + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1296.0   1303.8   -645.0   1290.0 

Coefficients: 
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.73622    2.87810 97.00000   0.256    0.799    
trtActive   -4.97544    0.99763 97.00000  -4.987 2.68e-06 ***
base         0.82712    0.09974 97.00000   8.293 6.37e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 34.4910 10.2957
6 10.2957 43.8530

Code

m1a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	trt	base	id	week
1	25.8	Placebo	31.9	1	4
2	23.8	Placebo	31.9	1	6
3	19.5	Active	26.5	2	4
4	21.0	Active	26.5	2	6
5	19.1	Active	25.8	3	4
6..199
200	12.3	Active	20.7	100	6

Code

summary(m1a)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	0.74	2.88	0.26	0.80
trtActive	−4.98	1.00	−4.99	<0.001
base	0.83	0.10	8.29	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model lead = trt base / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt / diff cl;
run;

Weekly TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0} + \beta_3time + \beta_4X\times time\]

Code

m1b <- mmrm(lead ~ trt + base + week + week:trt + 
               us(week|id),     
            data = tlc3_l %>% filter(week != '0'))

e1b <- emmeans(m1b, revpairwise ~ trt | week) %>% 
   as.data.frame()

summary(m1b)

mmrm fit

Formula:     lead ~ trt + base + week + week:trt + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1265.1   1272.9   -629.6   1259.1 

Coefficients: 
                 Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.91382    2.90529 100.46000   0.315    0.754    
trtActive        -7.35171    1.14443  97.99000  -6.424 4.80e-09 ***
base              0.82712    0.09974  97.00000   8.293 6.37e-13 ***
week6            -0.42400    0.94645  98.00000  -0.448    0.655    
trtActive:week6   5.67200    1.33848  98.00000   4.238 5.11e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 32.2160 13.4516
6 13.4516 39.4752

Code

m1b$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	trt	base	week	id
1	25.8	Placebo	31.9	4	1
2	23.8	Placebo	31.9	6	1
3	19.5	Active	26.5	4	2
4	21.0	Active	26.5	6	2
5	19.1	Active	25.8	4	3
6..199
200	12.3	Active	20.7	6	100

Code

summary(m1b)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	0.91	2.91	0.31	0.75
trtActive	−7.35	1.14	−6.42	<0.001
base	0.83	0.10	8.29	<0.001
week6	−0.42	0.95	−0.45	0.66
trtActive:week6	5.67	1.34	4.24	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model lead = trt base week week*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt*week / slicediff=week cl;
run;

Method 2: Repeated measures

Overall TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1X + \beta_2time + \beta_3time\times X\]

Code

m2a <- mmrm(lead ~ trt + follow + follow:trt + 
               us(week|id), 
            data = tlc3_l %>% mutate(follow = as.numeric(follow)-1)) 

e2a <- emmeans(m2a, revpairwise ~ trt | follow) %>% 
   as.data.frame()

summary(m2a)

mmrm fit

Formula:     lead ~ trt + follow + follow:trt + us(week | id)
Data:        %>%, tlc3_l, mutate(follow = as.numeric(follow) - 1) (used 300 
observations from 100 subjects with maximum 3 timepoints)
Covariance:  unstructured (6 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1900.0   1915.7   -944.0   1888.0 

Coefficients: 
                 Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)        28.019      0.705 98.000  39.743  < 2e-16 ***
trtActive          -1.763      0.997 98.000  -1.768   0.0801 .  
follow             -4.108      0.705 98.000  -5.827 7.21e-08 ***
trtActive:follow   -4.671      0.997 98.000  -4.685 9.03e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        0       4       6
0 25.0204 19.3844 22.5123
4 19.3844 49.1906 27.5832
6 22.5123 27.5832 63.7270

Code

m2a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	trt	follow	id	week
1	31.9	Placebo	0	1	0
2	25.8	Placebo	1	1	4
3	23.8	Placebo	1	1	6
4	26.5	Active	0	2	0
5	19.5	Active	1	2	4
6..299
300	12.3	Active	1	100	6

Code

summary(m2a)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	28.02	0.70	39.74	<0.001
trtActive	−1.76	1.00	−1.77	0.080
follow	−4.11	0.70	−5.83	<0.001
trtActive:follow	−4.67	1.00	−4.68	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 class id week (ref = first) follow (ref=first) trt (ref = last);
 model lead = trt follow follow*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt / diff cl;
run;

Weekly TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1X + \beta_2dummy\_time_1 + \beta_3dummy\_time_2 +\\ \beta_4dumm\_time_1\times X + \beta_5dummy\_time_2 \times X\]

Code

m2b <- mmrm(lead ~ trt + week + trt: week + us(week|id),
            data = tlc3_l)

e2b <- emmeans(m2b, revpairwise ~ trt | week) %>% 
   as.data.frame()

summary(m2b)

mmrm fit

Formula:     lead ~ trt + week + trt:week + us(week | id)
Data:        tlc3_l (used 300 observations from 100 subjects with maximum 3 
timepoints)
Covariance:  unstructured (6 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1869.1   1884.8   -928.6   1857.1 

Coefficients: 
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      27.9960     0.7069 98.0000  39.605  < 2e-16 ***
trtActive        -1.4560     0.9997 98.0000  -1.456    0.148    
week4            -3.9260     0.8132 98.0000  -4.828 5.09e-06 ***
week6            -4.3500     0.8887 98.0000  -4.895 3.87e-06 ***
trtActive:week4  -7.1000     1.1501 98.0000  -6.174 1.51e-08 ***
trtActive:week6  -1.4280     1.2567 98.0000  -1.136    0.259    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        0       4       6
0 24.9838 19.6480 22.0747
4 19.6480 47.3781 30.6207
6 22.0747 30.6207 58.6510

Code

m2b$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	trt	week	id
1	31.9	Placebo	0	1
2	25.8	Placebo	4	1
3	23.8	Placebo	6	1
4	26.5	Active	0	2
5	19.5	Active	4	2
6..299
300	12.3	Active	6	100

Code

summary(m2b)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	28.00	0.71	39.61	<0.001
trtActive	−1.46	1.00	−1.46	0.15
week4	−3.93	0.81	−4.83	<0.001
week6	−4.35	0.89	−4.90	<0.001
trtActive:week4	−7.10	1.15	−6.17	<0.001
trtActive:week6	−1.43	1.26	−1.14	0.26

Code

proc glimmix data=tlc3_l noclprint = 10;
 class id week (ref = first) follow (ref=first) trt (ref = last);
 model lead = trt follow follow*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt*follow / slicediff=follow cl;
run;

Method 2: Repeated measures wo/ treatment

Overall TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1time + \beta_2time\times X\]

Code

m2c <- mmrm(lead ~ follow + follow:trt + 
               us(week|id), 
            data = tlc3_l %>% mutate(follow = as.numeric(follow)-1))

e2f0 <- emmeans(m2c,            ~ trt | follow, at = list(follow = 0)) %>% as.data.frame()
e2f1 <- emmeans(m2c,revpairwise ~ trt | follow, at = list(follow = 1)) %>% as.data.frame()

e2c <- bind_rows(e2f0 %>% mutate(follow = '0'),
                 e2f1) 

summary(m2c)

mmrm fit

Formula:     lead ~ follow + follow:trt + us(week | id)
Data:        %>%, tlc3_l, mutate(follow = as.numeric(follow) - 1) (used 300 
observations from 100 subjects with maximum 3 timepoints)
Covariance:  unstructured (6 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1904.5   1920.1   -946.2   1892.5 

Coefficients: 
                 Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       27.2535     0.5026 99.0000  54.221  < 2e-16 ***
follow            -3.9753     0.6998 99.4800  -5.681 1.34e-07 ***
follow:trtActive  -4.9754     0.9820 98.0000  -5.067 1.91e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        0       4       6
0 25.2664 20.7527 21.1006
4 20.7527 51.2871 27.3801
6 21.1006 27.3801 61.2242

Code

m2c$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	follow	id	week	trt
1	31.9	0	1	0	Placebo
2	25.8	1	1	4	Placebo
3	23.8	1	1	6	Placebo
4	26.5	0	2	0	Active
5	19.5	1	2	4	Active
6..299
300	12.3	1	100	6	Active

Code

summary(m2c)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	27.25	0.50	54.22	<0.001
follow	−3.98	0.70	−5.68	<0.001
follow:trtActive	−4.98	0.98	−5.07	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 class id week (ref = first) trt (ref = last);
 model lead = follow follow*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
run;

Weekly TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} = \beta_0 + \beta_1dummytime_1 + \beta_2dummytime_2 + \\ \beta_3dummytime_1 \times X + \beta_4dummytime_2 \times X\]

Code

m2d <- mmrm(lead ~ I(week=='4') + I(week=='6') + 
               I(week=='4' & trt=='Active') + I(week=='6' & trt=='Active') + 
               us(week|id), 
            data = tlc3_l)

summary(m2d)

mmrm fit

Formula:     lead ~ I(week == "4") + I(week == "6") + I(week == "4" & trt ==  
    "Active") + I(week == "6" & trt == "Active") + us(week |      id)
Data:        tlc3_l (used 300 observations from 100 subjects with maximum 3 
timepoints)
Covariance:  unstructured (6 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1873.1   1888.7   -930.5   1861.1 

Coefficients: 
                                     Estimate Std. Error      df t value
(Intercept)                           27.2680     0.5027 99.0000  54.247
I(week == "4")TRUE                    -3.7705     0.8063 99.7200  -4.677
I(week == "6")TRUE                    -4.2652     0.8868 98.4200  -4.810
I(week == "4" & trt == "Active")TRUE  -7.4110     1.1301 98.0000  -6.558
I(week == "6" & trt == "Active")TRUE  -1.5975     1.2513 98.0000  -1.277
                                     Pr(>|t|)    
(Intercept)                           < 2e-16 ***
I(week == "4")TRUE                   9.17e-06 ***
I(week == "6")TRUE                   5.44e-06 ***
I(week == "4" & trt == "Active")TRUE 2.58e-09 ***
I(week == "6" & trt == "Active")TRUE    0.205    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        0       4       6
0 25.2668 19.8706 22.3248
4 19.8706 47.5531 30.8174
6 22.3248 30.8174 58.8718

Code

m2d$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	lead	I(week == "4")	I(week == "6")	I(week == "4" & trt == "Active")	I(week == "6" & trt == "Active")	id	week
1	31.9	FALSE	FALSE	FALSE	FALSE	1	0
2	25.8	TRUE	FALSE	FALSE	FALSE	1	4
3	23.8	FALSE	TRUE	FALSE	FALSE	1	6
4	26.5	FALSE	FALSE	FALSE	FALSE	2	0
5	19.5	TRUE	FALSE	TRUE	FALSE	2	4
6..299
300	12.3	FALSE	TRUE	FALSE	TRUE	100	6

Code

summary(m2d)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	27.27	0.50	54.25	<0.001
I(week == "4")TRUE	−3.77	0.81	−4.68	<0.001
I(week == "6")TRUE	−4.27	0.89	−4.81	<0.001
I(week == "4" & trt == "Active")TRUE	−7.41	1.13	−6.56	<0.001
I(week == "6" & trt == "Active")TRUE	−1.60	1.25	−1.28	0.20

Code

data tlc3_l;
 set tlc3_l;
 if week = "4" then week_4 = 1; else week_4 = 0;
 if week = "6" then week_6 = 1; else week_6 = 0;
run;

proc glimmix data=tlc3_l noclprint = 10;
 class id week (ref = first)  trt (ref = last);
 model lead = week_4 week_6 week_4*trt week_6*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
run;

Method 3: Analysis of changes (not adjusted)

Overall TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X\]

Code

m3a <- mmrm(diff ~ trt + 
               us(week|id), 
            data = tlc3_l %>% filter(week != '0'))

e3a <- emmeans(m3a, revpairwise ~ trt) %>% 
   as.data.frame()

summary(m3a)

mmrm fit

Formula:     diff ~ trt + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1296.2   1304.0   -645.1   1290.2 

Coefficients: 
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   -4.108      0.705 98.000  -5.827 7.21e-08 ***
trtActive     -4.671      0.997 98.000  -4.685 9.03e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 35.4435 10.7073
6 10.7073 43.7238

Code

m3a$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	diff	trt	id	week
1	-6.1	Placebo	1	4
2	-8.1	Placebo	1	6
3	-7.0	Active	2	4
4	-5.5	Active	2	6
5	-6.7	Active	3	4
6..199
200	-8.4	Active	100	6

Code

summary(m3a)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	−4.11	0.70	−5.83	<0.001
trtActive	−4.67	1.00	−4.68	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model diff = trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt / diff cl;
run;

Weekly TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X + \beta_2 time + \beta_3X \times time\]

Code

m3b <- mmrm(diff ~ trt + week + trt:week + 
               us(week|id), 
            data = tlc3_l %>% filter(week != '0'))

e3b <- emmeans(m3b, revpairwise ~ trt | week) %>% 
   as.data.frame()

summary(m3b)

mmrm fit

Formula:     diff ~ trt + week + trt:week + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1265.3   1273.1   -629.7   1259.3 

Coefficients: 
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      -3.9260     0.8132 98.0000  -4.828 5.09e-06 ***
trtActive        -7.1000     1.1501 98.0000  -6.174 1.51e-08 ***
week6            -0.4240     0.9464 98.0000  -0.448    0.655    
trtActive:week6   5.6720     1.3385 98.0000   4.238 5.11e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 33.0656 13.8818
6 13.8818 39.4858

Code

m3b$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	diff	trt	week	id
1	-6.1	Placebo	4	1
2	-8.1	Placebo	6	1
3	-7.0	Active	4	2
4	-5.5	Active	6	2
5	-6.7	Active	4	3
6..199
200	-8.4	Active	6	100

Code

summary(m3b)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	−3.93	0.81	−4.83	<0.001
trtActive	−7.10	1.15	−6.17	<0.001
week6	−0.42	0.95	−0.45	0.66
trtActive:week6	5.67	1.34	4.24	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model diff = trt week week*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt*week / slicediff=week cl;
run;

Method 3: Analysis of changes (adjusted)

Overall TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X +\beta_2Y_{t0}\]

Code

m3c <- mmrm(diff ~ trt + base + 
               us(week|id), 
            data = tlc3_l %>% filter(week != '0'))

e3c <- emmeans(m3c, revpairwise ~ trt) %>% 
   as.data.frame()

summary(m3c)

mmrm fit

Formula:     diff ~ trt + base + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1296.0   1303.8   -645.0   1290.0 

Coefficients: 
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.73622    2.87810 97.00000   0.256   0.7986    
trtActive   -4.97544    0.99763 97.00000  -4.987 2.68e-06 ***
base        -0.17288    0.09974 97.00000  -1.733   0.0862 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 34.4910 10.2957
6 10.2957 43.8530

Code

m3c$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	diff	trt	base	id	week
1	-6.1	Placebo	31.9	1	4
2	-8.1	Placebo	31.9	1	6
3	-7.0	Active	26.5	2	4
4	-5.5	Active	26.5	2	6
5	-6.7	Active	25.8	3	4
6..199
200	-8.4	Active	20.7	100	6

Code

summary(m3c)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	0.74	2.88	0.26	0.80
trtActive	−4.98	1.00	−4.99	<0.001
base	−0.17	0.10	−1.73	0.086

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model diff = trt base / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt / diff cl;
run;

Weekly TRT effect

EQ
MODEL
DATA
COEFF
SAS

\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X + \beta_2 Y_{t0} + \beta_3X \times time\]

Code

m3d <- mmrm(diff ~ trt + base + week + trt:week + 
               us(week|id), 
            data = tlc3_l %>% filter(week != '0'))

e3d <- emmeans(m3d, revpairwise ~ trt | week) %>% 
   as.data.frame()

summary(m3d)

mmrm fit

Formula:     diff ~ trt + base + week + trt:week + us(week | id)
Data:        %>%, tlc3_l, filter(week != "0") (used 200 observations from 100 
subjects with maximum 2 timepoints)
Covariance:  unstructured (3 variance parameters)
Method:      Satterthwaite
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  1265.1   1272.9   -629.6   1259.1 

Coefficients: 
                 Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       0.91382    2.90529 100.46000   0.315   0.7538    
trtActive        -7.35171    1.14443  97.99000  -6.424 4.80e-09 ***
base             -0.17288    0.09974  97.00000  -1.733   0.0862 .  
week6            -0.42400    0.94645  98.00000  -0.448   0.6551    
trtActive:week6   5.67200    1.33848  98.00000   4.238 5.11e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
        4       6
4 32.2160 13.4516
6 13.4516 39.4752

Code

m3d$tmb_data$full_frame %>%
   select(-`(weights)`) %>% 
   gt_preview(incl_rownums = TRUE)

	diff	trt	base	week	id
1	-6.1	Placebo	31.9	4	1
2	-8.1	Placebo	31.9	6	1
3	-7.0	Active	26.5	4	2
4	-5.5	Active	26.5	6	2
5	-6.7	Active	25.8	4	3
6..199
200	-8.4	Active	20.7	6	100

Code

summary(m3d)$coefficient %>% 
   as.data.frame() %>% 
   rownames_to_column() %>% 
   select(-df) %>% 
   gt(rowname_col = "rownames") %>% 
   fmt_number(columns = 2:4, decimals = 2) %>% 
   fmt(columns = 5, 
       fns = function(x) gtsummary::style_pvalue(x, digits = 2)) %>% 
   cols_align_decimal() %>% 
   cols_label(rowname = '')

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	0.91	2.91	0.31	0.75
trtActive	−7.35	1.14	−6.42	<0.001
base	−0.17	0.10	−1.73	0.086
week6	−0.42	0.95	−0.45	0.66
trtActive:week6	5.67	1.34	4.24	<0.001

Code

proc glimmix data=tlc3_l noclprint = 10;
 where week ^= 0;
 class id week (ref = first) trt (ref = last);
 model diff = trt base week week*trt / solution ddfm=satterthwaite;
 random _residual_ / subject=id type=un;
 lsmeans trt*week / slicediff=week cl;
run;

Overall TRT Effect (Combined Results)

	week	emmean	lower.CL	upper.CL	t.ratio	p.value
(1a) Longitudinal analysis of covariance	4/6	−5.0 ± 1.00	−7.0	−3.0	−4.99	<0.001
(2a) Repeated measures analysis	4/6	−6.4 ± 1.28	−9.0	−3.9	−5.02	<0.001
(2c) Repeated measures wo/ treatment	4/6	−5.0 ± 0.98	−6.9	−3.0	−5.07	<0.001
(3a) Analysis of changes (not adjusted)	4/6	−4.7 ± 1.00	−6.6	−2.7	−4.68	<0.001
(3c) Analysis of changes (adjusted)	4/6	−5.0 ± 1.00	−7.0	−3.0	−4.99	<0.001

(1a) overall TRT effect of -5.0 (\(\mu\) g/dL)
(2a) wo/ adjustment for baseline differences overestimation
(2b) w/ adjustment for baseline differences similar results to (1a)
(3a) wo/ adjustments for baseline difference underestimation
(3c) w/ adjustment for baseline differences exact same results

Mathematical equivalence between models

`Longitudinal analysis of covariance`	`Analysis of changes (adjusted)`
\[Y_{t} = \beta_0 + \beta_1X + \beta_2Y_{t0}\]	\[Y_{t} - Y_{t0} = \beta_0 + \beta_1X +\beta_2Y_{t0}\]
	\[Y_{t} = \beta_0 + \beta_1X +\beta_2Y_{t0} + Y_{t0}\]
	\[Y_{t} = \beta_0 + \beta_1X + (1 + \beta_2) Y_{t0}\]
\[Y_{t} = 0.74 - 4.98 \times X + 0.83 \times Y_{t0}\]	\[Y_{t} = 0.74 - 4.98 \times X - 0.17 \times Y_{t0}\]

TRT Effect (Combined Results Interactive)

TRT Effect w/ Missing (Combined Results Interactive)

Different ways to estimate treatment effects in randomised controlled trials

Generalized additive mixed models (GAMMs)

Generalized additive models (GAM)

GAM is a powerful and yet simple technique

Easy to interpret (Similar to GLM)
Flexible predictor functions can uncover hidden patterns in the data
Regularization of predictor functions helps avoid overfitting

(Generalized) Linear Models

\[y_i \sim \mathcal{N}(\mu_i, \sigma^2)\]

\[\mu_i = \beta_0 + \beta_1 \mathtt{month}_{i} + \beta_2 \mathtt{month}^2_{1i} + \cdots + \beta_j \mathtt{month}^j_{i}\]

Assumptions

linear effects of covariates are good approximation of the true effects

An additive model address the first of these

How is a GAM different?

\[\begin{align*} y_i &\sim \mathcal{D}(\mu_i, \theta) \\ \mathbb{E}(y_i) &= \mu_i = g(\eta_i)^{-1} \end{align*}\]

We model the mean of data as a sum of linear terms:

\[\eta_i = \beta_0 +\sum_j \color{red}{ \beta_j x_{ji}}\]

A GAM is a sum of smooth functions or smooths

\[\eta_i = \beta_0 + \sum_j \color{red}{f_j(x_{ji})}\]

Fitting a GAM in R

m3a <- gamm(anyexac2 ~ s(month, bs = 'cc'), 
              data = nejm, verbosePQL = FALSE,
              random  =  list(studyid = ~1), 
              family = binomial, 
              control = list(maxIter = 100), 
              method = 'REML')

m3  <- gamm(anyexac2 ~ group + s(month, bs = 'cc', by = group), 
             data = nejm, verbosePQL = FALSE,
             random  =  list(studyid = ~1), f
             amily = binomial, 
             control = list(maxIter = 100), 
             method = 'REML')

s() terms are smooths of one or more variables

bs = 'cc' terms specifies a cyclic cubic regression spline

\[\eta_i = \beta_0 + f(\mathtt{month}_i)\]

Seasonal Exacerbations Randomized double blind placebo control (RDBPC)

Month	N	Overall	Treatment		Difference¹	95% CI^1,2	p-value¹
Month	N	Overall	Control	Intervention	Difference¹	95% CI^1,2	p-value¹
Jan	377	9.3%	10.8%	7.9%	2.9%	-3.5%, 9.3%	0.43
Feb	380	10.5%	14.8%	6.3%	8.5%	1.9%, 15%	0.011
Mar	382	12.8%	14.2%	11.5%	2.8%	-4.5%, 10%	0.51
Apr	385	9.6%	14.1%	5.2%	9.0%	2.6%, 15%	0.005
May	380	8.9%	11.2%	6.7%	4.5%	-1.8%, 11%	0.18
Jun	375	5.1%	7.0%	3.2%	3.9%	-1.1%, 8.8%	0.14
Jul	372	4.0%	5.5%	2.6%	2.8%	-1.7%, 7.4%	0.26
Aug	375	6.1%	5.9%	6.3%	-0.44%	-5.7%, 4.9%	>0.99
Sep	370	9.5%	13.6%	5.4%	8.2%	1.8%, 15%	0.012
Oct	369	11.1%	13.7%	8.6%	5.1%	-1.9%, 12%	0.17
Nov	371	8.6%	13.0%	4.3%	8.7%	2.5%, 15%	0.005
Dec	374	9.9%	10.9%	8.9%	1.9%	-4.7%, 8.5%	0.65
¹ Two sample test for equality of proportions
² CI = Confidence Interval

Additive Model

Code

m3a <- gamm(anyexac2 ~ s(month, bs = 'cc'), data = nejm, verbosePQL = FALSE,
            random  =  list(studyid = ~1), family = binomial, control = list(maxIter = 100), method = 'REML')

m3  <- gamm(anyexac2 ~ group + s(month, bs = 'cc', by = group), data = nejm, verbosePQL = FALSE,
            random  =  list(studyid = ~1), family = binomial, control = list(maxIter = 100), method = 'REML')

OVERALL

BY TREATMENT

Derivatives

Code

# DERIVATVES
f3a_d1 <- derivatives(m3a,
                      level = 0.68,
                      type = 'central',
                      interval = "simultaneous",
                      order = 1) %>% 
   draw() +
   geom_pointless(aes(color = after_stat(location)),
                  location = c("min","max"),
                  size = 3,
                  color = 'gray25') +
   geom_hline(yintercept = 0, color = 'gray50') + 
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   scale_y_continuous(limits = c(-0.8, 0.8) ) +
   labs(x  =  "Month",
        y  = "1st Derivative (Rate of Change)",
        title = NULL) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

f3a_d2 <- derivatives(m3a,
                      level = 0.68,
                      type = 'central',
                      interval = "simultaneous",
                      order = 2) %>% 
   draw() +
   geom_pointless(aes(color = after_stat(location)),
                  location = c("min","max"),
                  size = 3,
                  color = 'gray25') +
   geom_hline(yintercept = 0, color = 'gray50') + 
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   scale_y_continuous(limits = c(-1.2, 1.2) ) +
   labs(x  =  "Month",
        y  = "2nd Derivative (Inflection Point)",
        title = NULL) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

# DERIVATVES
d3d1 <- derivatives(m3,
                    level = 0.68,
                    type = 'central',
                    interval = "simultaneous",
                    order = 1) %>% 
   mutate(group = factor(smooth, labels = c("Control","Intervention")))

f3_d1 <- ggplot(data  =  d3d1,
                aes(x  =  data, y  =  derivative, color  =  group, label = group)) +
   geom_pointless(aes(color = after_stat(location)),
                  location = c("min","max"),
                  color = 'gray25',
                  size = 3) + 
   facet_wrap(~ group) +
   geom_line(linewidth  =  1) +
   geom_ribbon(aes( ymin  =  lower, ymax  =  upper, 
                    color  =  NULL, fill  =  smooth, alpha  =  0.5))+
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   scale_y_continuous(limits = c(-0.8, 0.8) ) +
   scale_color_brewer(type  =  'qual', palette  =  'Dark2') + 
   scale_fill_brewer(type  =  'qual', palette  =  'Dark2') + 
   labs(x  =  "Month",
        y  =  "1st Derivative (Rate of Change",
        title  =  ) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

d3d2 <- derivatives(m3,
                    level = 0.68,
                    type = 'central',
                    interval = "simultaneous",
                    order = 2) %>% 
   mutate(group = factor(smooth, labels = c("Control","Intervention")))

f3_d2 <- ggplot(data  =  d3d2,
                aes(x  =  data, y  =  derivative, color  =  group, label = group)) +
   geom_pointless(aes(color = after_stat(location)),
                  location = c("min","max"),
                  color = 'gray25',
                  size = 3) + 
   facet_wrap(~ group) +
   geom_line(linewidth  =  1) +
   geom_ribbon(aes( ymin  =  lower, ymax  =  upper, 
                    color  =  NULL, fill  =  smooth, alpha  =  0.5))+
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   scale_y_continuous(limits = c(-1.2, 1.2) ) +
   scale_color_brewer(type  =  'qual', palette  =  'Dark2') + 
   scale_fill_brewer(type  =  'qual', palette  =  'Dark2') + 
   labs(x  =  "Month",
        y  =  "2nd Derivative (Inflection Point)",
        title  =  ) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

design <- "
144
255
366
"

f3a + 
   f3a_d1 + 
   f3a_d2 + 
   (f3+facet_wrap(~group) + theme(strip.background = element_blank(),strip.text.x = element_blank())) + 
   (f3_d1 + theme(strip.background = element_blank(),strip.text.x = element_blank())) + 
   (f3_d2 + theme(strip.background = element_blank(),strip.text.x = element_blank())) +
   plot_layout(design = design)

Smooth Differences

Code

d2s1 <- smooth_samples(m3$gam,
                       n = 500,
                       data = d1,
                       seed = 123) %>% 
   rename(month = .x1) %>% 
   group_by(draw, month) %>% 
   summarise(diff = value[1]-value[2]) %>% 
   mutate(diff = (diff + 0.35)/10 )


f3_s1 <- ggplot(data  =  d2s1,
                aes(x  =  month, y  =  diff, group = draw, label = 'Difference')) +
   geom_line(linewidth  =  0.25, alpha = 0.25) +
   geom_hline(yintercept = 0, color = 'gray50') +
   scale_y_continuous(limits  =  c(-0.05, 0.12),
                      minor_breaks  =  FALSE,
                      expand  =  c(0, 0.01)) + 
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   labs(x  =  "Month",
        y  =  "Difference between smooths",
        title  =  ) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

s1b <- difference_smooths(m3$gam, 
                          smooth = "s(month)",
                          ci_level = 0.95) %>% 
   mutate(across(c(diff, lower, upper), ~(.x + 0.35)/10)) 

f3_s2 <- ggplot(data  =  s1b,
                aes(x  =  month, y  =  diff,  color = (lower>0) ) ) +
   geom_hline(yintercept = 0, color = 'gray50') +
   geom_line(aes(group = 1), linewidth  = 2) +
   geom_ribbon(data = s1b,
               inherit.aes = FALSE,
               aes(x = month, ymin  = lower, ymax  =  upper,  alpha  =  0.5))+
   scale_y_continuous(limits  =  c(-0.05, 0.12),
                      minor_breaks  =  FALSE,
                      expand  =  c(0, 0.01)) + 
   scale_x_continuous(breaks  =  seq(1, 12, 1),
                      labels  =  c("J","F","M","A","M","J","J","A","S","O","N","D"),
                      expand  =  c(0.01, 0.01),
                      minor_breaks  =  FALSE) +
   labs(x  =  "Month",
        y  =  "Difference between smooths",
        title  =  NULL) +
   theme_bw(base_size = 20) +
   theme(legend.position  =  "none",
         panel.grid.minor  =  element_blank())

f3/f3_s1/f3_s2

Data Visualizations for Clinical Trials Reporting

Visualizations for the Special Interest Group at Statisticians in the Pharmaceutical Industry (PSI)

All Analyses @agstn/WW

Thanks!

R Session Information for reproducibility

Setting	Value
R environment
version	R version 4.2.0 (2022-04-22 ucrt)
os	Windows 10 x64 (build 17763)
system	x86_64, mingw32
ui	RTerm
language	(EN)
collate	English_United States.1252
ctype	English_United States.1252
tz	America/Los_Angeles
date	2023-01-31
pandoc	2.19.2 @ C:/R/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

Package	Loaded version	Date
R Packages
broom	1.0.3	2023-01-25
dplyr	1.0.10	2022-09-01
emmeans	1.8.4-1	2023-01-17
forcats	0.5.2	2022-08-19
geomtextpath	0.1.1	2022-08-30
ggbeeswarm	0.7.1	2022-12-16
ggplot2	3.4.0	2022-11-04
ggpointless	0.0.3	2022-08-25
gratia	0.7.3	2022-05-09
gt	0.8.0	2022-11-16
gtreg	0.2.0	2022-10-18
gtsummary	1.7.0	2023-01-13
labelled	2.10.0	2022-09-14
MASS	7.3-58.2	2023-01-23
mgcv	1.8-41	2022-10-21
mmrm	0.2.2	2022-12-20
multcomp	1.4-20	2022-08-07
mvtnorm	1.1-3	2021-10-08
nlme	3.1-161	2022-12-15
patchwork	1.1.2	2022-08-19
purrr	1.0.1	2023-01-10
reactable	0.4.3	2023-01-07
reactablefmtr	2.1.0	2022-06-27
readr	2.1.3	2022-10-01
rio	0.5.29	2021-11-22
scatterPlotMatrix	0.2.0	2022-04-11
stringr	1.5.0	2022-12-02
survival	3.5-0	2023-01-09
TH.data	1.1-1	2022-04-26
tibble	3.1.8	2022-07-22
tidyr	1.3.0	2023-01-24
tidyverse	1.3.2	2022-07-18
trelliscopejs	0.2.6	2021-02-01

Estimate Treatment Effects RCT

Outline

How to estimate treatment effects in longitudinal RCT’s

Problem

Data

Models R &SAS

Results

How to extend these estimation

How to aggregate & communicate these results

Different ways to estimate treatment effects in randomised controlled trials

Example dataset

Table (wo/missing)

Tables (w/missing)

Figures

Method 1: Longitudinal analysis of covariance

Overall TRT effect

Weekly TRT effect

Method 2: Repeated measures

Overall TRT effect

Weekly TRT effect

Method 2: Repeated measures wo/ treatment

Overall TRT effect

Weekly TRT effect

Method 3: Analysis of changes (not adjusted)

Overall TRT effect

Weekly TRT effect

Method 3: Analysis of changes (adjusted)

Overall TRT effect

Weekly TRT effect

Overall TRT Effect (Combined Results)

Mathematical equivalence between models

TRT Effect (Combined Results Interactive)

TRT Effect w/ Missing (Combined Results Interactive)

Generalized additive mixed models (GAMMs)

Generalized additive models (GAM)

(Generalized) Linear Models

How is a GAM different?

Seasonal Exacerbations Randomized double blind placebo control (RDBPC)

Additive Model

OVERALL

BY TREATMENT

Derivatives

Smooth Differences

Data Visualizations for Clinical Trials Reporting

Thanks!

R Session Information for reproducibility